% T2_FSE: Estimates parameters and errorbars for T2 relaxation 
% Signal model = A*exp(-TE/T2)+C

% Execute by hitting F5

% After selection of the image it will pop up and two ROIs need to be drawn: First a ROI with only background noise and 
% second the ROI where the mapping is desired.

% Output are dicom images of the Amap, R2map, R2errormap, T1map and T1errormap and Cmap, both for Maximum Likelihood
% and Least Squares estimation. They are saved in the directory of the first image of the serie. 

% IMPORTANT: In the maps the values of the R2 and R2 errormap are multiplied with 10000 to get into sensible gray values.

% The errors in the errormap are the lower bounds on the standard deviations calculated with
% the Cramer-Rao lower bounds.

% Created by Henk Smit, EMC, 01-2011 based on the work by Dirk Poot,
% University of Antwerp, 13-8-2007. 

clear all

% switch for voxel/roi/roimean type of estimation
% known bugs: 1 (all voxels) only succeeds when every TE only occurs once
% the LS result with 1 (all voxels) is not accurate
% the std with 2(roi mean) is too high. Use 1 (all voxels) instead.
roibased = 2; % 0=voxel, 1=roi all voxels, 2=roi mean

% choose if you want a T2 or a R2 map overlay image
overlaytype = 0; % 0=T2, 1=R2

% with matlabvs > 7 a dialog window can be used 
% the roibased = 0/1 option
% roivoxel = VoxelRoi; 
% roibased = strcmp('ROI',roivoxel);

% fill after how many images the same cardiac phase is reached
repeatphase = 1;

% which should be the first point to take into account
fitfrom = 1; 

% output values as a txt file
txtfile = true;

[file,chosendir]=uigetfile('*.dcm');
imagenr = str2num(file(1:end-4));
cd(chosendir);
d=dir;
nfol=(length(d)-2)/repeatphase;
cd ..
studydir = pwd; %save maps in dir of chosen image

%read in the image that is picked by the user
chosenimage = dicomread(dicominfo(fullfile(chosendir,file)));

%get the dicom info
for i=1:nfol   
    %henk orig fitdata(i).folder=char(d((i-1)*repeatphase+2+imagenr).name); 
    fitdata(i).folder=char(d((i-1)*repeatphase+2+1).name); 
    info = dicominfo(fullfile(chosendir,fitdata(i).folder));
    fitdata(i).image = double(dicomread(info));
    fitdata(i).TE=info.EchoTime;
    fitdata(i).TR=info.RepetitionTime+99999; %TR is not correctly in header, should look it up. It's long compared to T1 anyway, so almost neglectible
    fitdata(i).Angle=info.FlipAngle;
    fitdata(i).info=info;
end


% % file logistics
% [file,chosendir]=uigetfile('*.dcm');
% cd(chosendir);
% cd ..
% studydir=pwd;
% cd(studydir)
% d=dir;
% nfol=length(d);
% fitdata(1).folder=char(d(1+2).name);
% 
% %get the dicom info
% for i=1:nfol-2
%     fitdata(i).folder=char(d(i+2).name);
%     info = dicominfo([studydir,'\',fitdata(i).folder,'\',file]);
%     fitdata(i).image = double(dicomread(info));
%     fitdata(i).TE=info.EchoTime;
%     fitdata(i).TR=info.RepetitionTime;
%     fitdata(i).Angle=info.FlipAngle;
%     fitdata(i).info=info;
% end
% 
% %read in the image that is picked by the user
% chosenimage = dicomread(dicominfo(fullfile(chosendir,file)));

%determine the sigma of the noise
imshow(chosenimage,[]);
text(35,120,'Select a region in the background','Color','yellow')
M=roipoly;
close;

%read in data to fitdata struct
for i=1:nfol-2
    D=fitdata(i).image.*M;
    sdyd(i)=std(nonzeros(D));
    maxsd=max(sdyd);
end

clear D;

%conversion factor from rayleigh distribution
CRLBsigma=1.527*maxsd; 
% CRLBsigma=5;

%draw rois
% cd(studydir)
% d=dir;
% nfol=length(d);
for i=1:nfol
    if i==1
        imshow(chosenimage,[]);
        text(45,350,'Select the region of interest','Color','yellow')
% %when ROI is to be drawn here:
        M=roipoly;
        close;
%when ROI is imported
    %   M = maskimage ./ maskimage;
%do the whole image   
%       M=ones(256,256);
    
    end

    fitdata(i).image=fitdata(i).image.*M;
end


TE=zeros(nfol,1);

for k=1:nfol
    TE(k)=double(fitdata(k).TE);
end

%initialize&allocate
TEprime=TE;
yd=zeros(size(TE,1),size(TE,2));

if ~roibased
    [row,column]=find(fitdata(2).image);
    nrvoxels=size(row,1);
end;


dims=size(fitdata(1).image);
R2=zeros(dims(1),dims(2));

dims=size(fitdata(1).image);
T2=zeros(dims(1),dims(2));
R2=T2;
A=T2;
T2ls=T2;
C=T2;
R2ls=T2;
Als=T2;
Cls=T2;
STDCRT2=T2;
STDCRR2=T2;
STDCRA=T2;


% setting fit and ML options
opt = optimset('fminunc');
opt = optimset(opt,'LargeScale','off','gradObj','on','DerivativeCheck','off','Display','off','MaxIter',40,'Hessian','off','TolFun',1e-11,'TolX',1e-9,'MaxFunEvals',1000);
s = fitoptions('Method','NonLinearLeastSquares','Lower',[-10000,-10000,-100],'Upper',[50000,10000,10000],'Startpoint',[1000,0,0],'Maxiter',80,'TolFun',10^-9,'TolX',10^-9,'Display','off');
f = fittype('a*exp(-x*b)+c','options',s);
fitrange = 1:size(TE);
% for FMINCON opt =
% optimset(opt,'gradObj','on','Display','off','MaxIter',40,'Hessian','off','TolFun',1e-8,'TolX',1e-7);

% logistics for roibased all voxels
 if roibased == 1
     yd=double(fitdata(1).image);
     for k=2:nfol
            yd=[yd double(fitdata(k).image)];
     end
     yd=reshape(yd,size(yd,1)*size(yd,2),1);
     TE=repmat(TE,1,size(yd,1)/(nfol));
     TE=TE';
     TE=reshape(TE,size(yd,1),1);
     [nonzerox,nonzeroy]=find(yd);
     yd=yd(nonzerox,1);
     TE=TE(nonzerox,1);
     nrvoxels=1;
     
% logistics for roibased mean of voxels
 elseif roibased ==2
    yd=mean(double(nonzeros(fitdata(1).image)));
         for k=2:nfol
            yd=[yd mean(double(nonzeros(fitdata(k).image)))];
         end  
        yd=yd';
        nrvoxels = 1;      
 end
    
% actual estimation
for i=1:nrvoxels
    disp(['Calculating pixel: ' num2str(i) '/' num2str(nrvoxels)]);
    if ~roibased   
        if fitdata(1).image(row(i),column(i)) > 0;
                for k=1:nfol
                    yd(k)=double(fitdata(k).image(row(i),column(i)));
                end

                [LS,ML]=T2ABCComputePar(yd,TE,0,CRLBsigma,opt,s,f,fitrange);
                [CR]=T2ABCCramerRao(ML,TE,CRLBsigma);

                A(row(i),column(i)) = ML(1,1);
                R2(row(i),column(i)) = ML(2,1);
                T2(row(i),column(i)) = 1/ML(2,1);
                C(row(i),column(i)) = ML(3,1);
                
                Als(row(i),column(i)) = LS(1,1);
                R2ls(row(i),column(i)) = LS(2,1);
                T2ls(row(i),column(i)) = 1/LS(2,1);
                Cls(row(i),column(i)) = LS(3,1);

                if CR(1,1)>0 && ~isnan(CR(1,1))
                    STDCRA(row(i),column(i)) = sqrt(CR(1,1));
                else
                    STDCRA(row(i),column(i)) = 10000;
                end
                if CR(2,2)>0 && ~isnan(CR(2,2))
                    STDCRR2(row(i),column(i)) = sqrt(CR(2,2));
                    STDCRT2(row(i),column(i)) = sqrt(CR(2,2))/(ML(2,1)^2);
                else
                    STDCRR2(row(i),column(i)) = 5;
                    STDCRT2(row(i),column(i)) = 10000;
                end

            else
                A(row(i),column(i))= 0;
                T2(row(i),column(i)) = 0;
                R2(row(i),column(i)) = 0;
                Als(row(i),column(i))= 0;
                T2ls(row(i),column(i)) = 0;
                R2ls(row(i),column(i)) = 0;
                STDCRR2(row(i),column(i)) = 0;
                STDCRT2(row(i),column(i)) = 0;
        end  
    else
        [LS,ML]=T2ABCComputePar(yd,TE,0,CRLBsigma,opt,s,f,fitrange);
        [CR]=T2ABCCramerRao(ML,TE,CRLBsigma);
        A = ML(1,1);
        T2 = 1/ML(2,1);
        R2 = ML(2,1);
        C = ML(3,1);
        Als = LS(1,1);
        T2ls = 1/LS(2,1);
        R2ls = LS(2,1);
        Cls = LS(3,1);
        if CR(1,1)>0 && ~isnan(CR(1,1))
            STDCRA = sqrt(CR(1,1));
        else
            STDCRA = 10000;
        end
        if(CR(2,2)>0 && ~isnan(CR(2,2)))
            STDCRR1 = sqrt(CR(2,2));
            STDCRT1 = sqrt(CR(2,2))/(ML(2,1)^2); %error propagation dT/T = dR/R
        else
            STDCRR1 = 5;
            STDCRT1 = 10000;
        end
    end
end

%write result in dcm files. values*10000 to display the values in
%milliseconds, sqrt of CRR2 is to go to stds instead of vars.

if ~roibased
    dicomwrite(int16(10000*R2),[studydir,'\',fitdata(1).folder,'\',file(1:end-4),'R2map','.dcm'],fitdata(1).info);
    dicomwrite(int16(T2),[studydir,'\',fitdata(1).folder,'\',file(1:end-4),'T2map','.dcm'],fitdata(1).info);
    dicomwrite(int16(A),[studydir,'\',fitdata(1).folder,'\',file(1:end-4),'Amap','.dcm'],fitdata(1).info);
    dicomwrite(int16(C),[studydir,'\',fitdata(1).folder,'\',file(1:end-4),'Cmap','.dcm'],fitdata(1).info);
    dicomwrite(int16(10000*STDCRR2),[studydir,'\',fitdata(1).folder,'\',file(1:end-4),'R2errormap','.dcm'],fitdata(1).info);
    dicomwrite(int16(STDCRT2),[studydir,'\',fitdata(1).folder,'\',file(1:end-4),'T2errormap','.dcm'],fitdata(1).info);
    
    dicomwrite(int16(10000*R2ls),[studydir,'\',fitdata(1).folder,'\',file(1:end-4),'R2mapLS','.dcm'],fitdata(1).info);
    dicomwrite(int16(T2ls),[studydir,'\',fitdata(1).folder,'\',file(1:end-4),'T2mapLS','.dcm'],fitdata(1).info);
    dicomwrite(int16(Als),[studydir,'\',fitdata(1).folder,'\',file(1:end-4),'AmapLS','.dcm'],fitdata(1).info);
    dicomwrite(int16(Cls),[studydir,'\',fitdata(1).folder,'\',file(1:end-4),'CmapLS','.dcm'],fitdata(1).info);

    
    disp(['Resulting maps are stored in ' fullfile(studydir,fitdata(1).folder) ])
    
    if ~overlaytype
        %T2 overlay image. The "[# #]" values are the range of the T2 map.
        %The range of the image is set as [0 mean+6*std]
        %The range of the map is set as [0 3*median]
        imagecolorrange = [0 mean(nonzeros(chosenimage))+6*std(single(nonzeros(chosenimage)))];
        mapcolorrange = [0 2.5*median(nonzeros(T2))];
        figure('Name','Overlay with T2 map (milliseconds)', 'position', [800 600 500 500]);
        sc(T2ls,[mapcolorrange], jet,sc(chosenimage, [imagecolorrange],gray),T2==0);
        colorbar('south','XColor','white');
    else
        %R2 overlay image. The "[# #]" values are the range of the R2 map.
        imagecolorrange = [0 mean(nonzeros(chosenimage))+6*std(single(nonzeros(chosenimage)))];
        mapcolorrange = [0 3*median(nonzeros(R2))];
        figure('Name','Overlay with R2 map (R2*10000, 1/milliseconds)', 'position', [800 600 500 500]);
        sc(R2,[0 max(nonzeros(R2))], jet,sc(chosenimage, [imagecolorrange],gray),R2==0);
        colorbar('south','XColor','white');
    end
    
        disp(['Median T2 (ms) = ' num2str(median(nonzeros(T2))) ' +- ' num2str(median(nonzeros(STDCRT2)))])
   
end

%plotting for roibased fit
if roibased
    scatter(TE,yd,'*k')
    hold on
    te=floor(min(TE)-5):1:round(max(TE)+15);
    yy=LS(1,1)*exp(-te*LS(2,1))+LS(3,1);
    yy2=ML(1,1)*exp(-te*ML(2,1))+ML(3,1);
%     yy2=353*exp(-te*0.12971)+79;
%     plot(te,yy,'b');
%     scatter(TE,yd2,'*r')
    plot(te,yy2,'--k');
    hold on
    
    disp(['Model = A * Exp(-R2*TE) + B'])
    disp(['Results: Least Squares estimate'])
    disp(['R2 = ' num2str(R2ls)])
    disp(['T2 = ' num2str(T2ls)])
    disp(['A = ' num2str(Als)])
    disp(['B = ' num2str(Cls)])
    disp(['Results: Maximum Likelihood estimate +- Cramer-Rao Lower Bound on standard deviation:'])
    disp(['R2 = ' num2str(ML(2,1)) ' +- ' num2str(sqrt(CR(2,2)))])
    disp(['T2 = ' num2str(1/ML(2,1)) ' +- ' num2str(sqrt(CR(2,2))/(ML(2,1)^2))])
    disp(['A = ' num2str(ML(1,1)) ' +- ' num2str(sqrt(CR(1,1)))])
    disp(['B = ' num2str(ML(3,1)) ' +- ' num2str(sqrt(CR(3,3)))])


end;

%lx=nonzeros(EigVec(1,1,:)).*sqrt(nonzeros(EigVal(1,1,:)));
%ux=nonzeros(EigVec(1,2,:)).*sqrt(nonzeros(EigVal(2,2,:)));
%ly=nonzeros(EigVec(2,1,:)).*sqrt(nonzeros(EigVal(1,1,:)));
%uy=nonzeros(EigVec(2,2,:)).*sqrt(nonzeros(EigVal(2,2,:)));

%To show an A vs. R2 map with errorbars(either parallel with x and y axis
%or in the direction of the eigenvectors of the covariancematrix) enable one of the following lines: 

%errorbarxy(nonzeros(A)',nonzeros(R2)',sqrt(nonzeros(CRA))',(sqrt(nonzeros(CRR2)))',[],[],'w')
%errorbarxyoblique(nonzeros(A)',nonzeros(R2)',lx',ly',ux',uy','w')
%hold on
%scatter(nonzeros(A),nonzeros(R2),'.k')
 

